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The precise determination of the universality classes in self-organized critical phenomena (SOC) 
is still an unsolved problem. Different SOC models like sandpile, linear interface depinning, and the 
Barkhausen effect have been investigated independently. In the present work we demonstrate that 
these models can all be mapped into a linear interface depinning model driven at constant rate. 
The model is found to belong to the universality class of constant force linear interface depinning 
above the depinning transition, with an upper critical dimension d c = 4. Results are compared with 
numerical simulations, experiments, and previous theoretical works reported in the literature. In 
this way we demonstrate a precise connection between different SOC models which display the same 
universal behavior. 
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I. INTRODUCTION 



About one decade ago Bak, Tang and Wiesenfeld 
(BTW) introduced the idea of self-organized critical- 
ity (SOC) pj] to describe the critical behavior of a 
vast class of driven dissipative systems. The dynam- 
ics of SOC systems is characterized by long periods 
of quasi-equilibrium alternating with sudden rearrange- 
ments (avalanches) which may expand all the system. 
Avalanche dynamics has been observed in granular sand- 
piles Jj],§, the Barkhausen effect superconducting 
vortex piles || , earthquakes || , etc. 

As usual the mean-field (MF) theory gave the first in- 
sight into SOC phenomena [^fj|-p^|. We also count with 
some exact results for Abelian sandpile models p3| . [T4| . 
real space |Te| ] and momentum space [ fl6| renormaliza- 
tion group (RG) calculations, and singular diffusion equa- 
tions |l7| . More recently the analogy with other non- 
equilibrium systems like systems with many absorbing 
states jl8],[l9| , directed percolation [p0[-p2[ , and linear in- 
terface depinning (LID) |2^-^rj] has been exploited. This 
accumulated experience allowed a more precise formula- 
tion of SOC phenomena. It is now well known that the 
absence of contro 1 parameters, as it was believed at the 
early state, is not completely true. Some hidden param- 
eters like the driving rate should be fine-tuned in order 
to reach the critical state. 

However the precise determination of the universality 
classes in SOC phenomena is still an unsolved problem. 
For instance, it is not clear if deterministic and stochas- 
tic sandpile models belong to the same universality class. 
Some numerical simulations p7|,j28f , real space RG calcu- 
lations [H| and field theory |19(| pointed to the affirma- 
tive answer but other numerical simulations f2fj[| disagree. 
The problem becomes more difficult since in determinis- 
tic models multifractal scaling, instead of finite size scal- 
ing, is satisfied |3(|. On the other hand, some authors 
have pointed out the existence of universal behavior be- 
tween sandpile and LID models |24|,^6|]. Paczuski and 
Boettcher Q mapped the random slope threshold sand- 



pile model |31|] into a model of LID where the interface is 
pulled at one end. More recently we have extended this 
analogy to critical height models pq| , which are mapped 
into a LID model driven at constant rate. In this case 
we have obtain ed the complete set of scaling exponents 
using previous RG analysis for constant force LID mod- 
els. Lauritsen and Alava |2(J have also considered the 
analogy between critical height models and LID models. 
However in their work the way in which driving and dis- 
sipation take place are not well specified. The analysis is 
limited to a qualitative level and no rigorous determina- 
tion of the scaling exponents was provided, thought they 
count with the previous work by us p5[ . 

The way in which the interface is driven seems to 
be an important condition to obtain an ordinary depin- 
ning transition or SOC behavior. The interface may be 
driven either by extremal dynamics [|32|-|34"| , by constant 
force |35|-|37j or by constant rate |25|. While constant 
force models has been extensively studied in the liter- 
ature extremal dynamics and constant rate models are 
less known. However, in the last years, extremal dynam- 
ics and constant velocity models have gained more atten- 
tion due to its relation with the theory of self-organized 
criticality (SOC). However, constant rate, and not con- 
stant force or constant velocity, is the predominant driv- 
ing mechanism in experiments where avalanche dynam- 
ics is observed. Examples are found in granular sandpiles 
where the number of grains added to the system increases 
linearly with time Q , in Barkhausen and superconduct- 
ing avalanches experiments where the applied magnetic 
field is increased at constant rate, and in earthquake dy- 
namics where some tectonic plate dynamics gives a con- 
stant rate of stress increase leading to a stick-slip motion 
of other plates. 

In the present work we investigate the existence of uni- 
versality in different SOC models proposed in the litera- 
ture, which can be mapped into a LID model driven at 
constant rate. The LID model driven at constant rate is 
taken here as the prototype of this class of SOC models 
because it allows some analytical treatment, which is very 
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difficult for instance in sandpile models. The existence of 
a different driving mechanism does not change the univer- 
sality class; we show that the constant rate model belongs 
to the universality class of constant force LID model, with 
an upper critical dimension d c — 4. Contrary to previous 
works, our analysis is not limited to a qualitative level 
but it is supported by a complete mapping of different 
SOC models into a LID model. The complete set of scal- 
ing exponents, including the avalanche exponents, are 
computed using previous RG analysis developed for con- 
stant force LID models. Part of these results are already 
contained in [p5[ . 

The paper is organized as follows. In section |H| wc 
introduce the LID model driven at constant force and 
review some of the results obtained from RG analysis. 
The LID model driven at constant force has been already 
investigated in the literature, we just present some im- 
portant aspects which are determinant for the analysis 
developed in the next section. In section III wc consider 
the constant rate LID model and illustrate how different 
SOC models can be mapped into this model. The phase 
diagram of this model is investigated in section [V. The 



comparison of our results with experiments, numerical 
simulations and other theoretical approaches proposed 
in the literature is discussed in section [v|. Finally the 



summary and conclusions are given in section VI 



II. CONSTANT FORCE LID MODELS 

There are many physical situations where one is in- 
terested in an interface dynamics. Examples are found 
in the domain wall dynamics [§,[§, the displacement of 
one fluid by another inside a porous media |58|,[59]], con- 
tact line depinning [Q, and more In general a d- 
dimensional self-afhne interface, described by a single- 
valued function h{x,t), evolves in a (d + l)-dimensional 
medium. Usually some kind of disorder affects the mo- 
tion of the interface leading to its roughening. Earlier 
studies [[l2| focuse on time-independent uncorrelated dis- 
order but most recent studies analyze the motion of inter- 
faces under quenched disorder |[9|^[l3|,[l4| . In the presence 
of quenched disorder and constant force driving two uni- 
versality classes has been found f45| . One is described 
by the Kardar-Parisi-Zhang equation (42| with quenched 
noise. In this case the interface is pinned by paths on 
a directed percolation cluster of pinning sites p(|. The 
second class is described by the Edwards- Wilkinson equa- 
tion p7j | with quenched noise, usually known as LID mod- 
els. 

In LID models the interface height satisfies the equa- 
tion of motion 



Xd t h = rV 2 /i + F + r)(x,h), 



(1) 



where dt denotes the partial time derivative and V is a 
rf-dimensional Laplacian. A is a viscosity coefficient, T 
the surface tension and F is a constant force acting on 



the interface. r](x, h) is a random pinning force associ- 
ated with the existence of random pinning centers in the 
(d+l)-dimensional environment. In general it is assumed 
that r)(x, h) is a Gaussian noise uncorrelated in the space 
x, with zero mean and noise correlator 



(r)(x, h)r)(x', h')) = 5 d (x - x')A(h - ti). 



(2) 



where A(— h) = A(h) is a symmetric function, with a 
fast decay to zero beyond the characteristic length of the 
pinning centers. 

The existence of pinning centers will carry as a conse- 
quence that the interface will be not smooth. The rough- 
ness of the interface is characterized by the height-height 
correlation function 



([ y (M)-y(o,o)n~M 2C <?(i/l* 



(3) 



where £ and z are the roughness and dynamic scal- 
ing exponents, respectively, and g{x) a scaling function 
with the asymptotic behaviors g(x) ~ 1 for x <C 1 and 
g(x) ~ a; 2 ^/ z for x 3> 1. The scaling exponents z and 
£ are related with the scaling behavior of the interface 
fluctuations just at the critical state and therefore we ex- 
pecct that they are independent on the way the system 
is driven out from the critical state. 

A depinning transition takes place at certain critical 
force F c determined by the disorder. For F < F c the in- 
terface is pinned after certain finite time while for F > F c 
it moves with finite average velocity v. Above F c one 
looks for a solution of eq. (|l|) as an expansion around 
the flat co-moving interface vt, i.e. h(x, t) = vt + y(x, t), 
resulting the equation of motion for y(x,t) 



A<9 t y 



rV 2 j/ + F — Aw + r)(x, vt + y). 



(4) 



v is obtained self-consistently using the constraint 
(y(x, t)) = 0. The noise term in eq. (|J) makes this equa- 
tion non-linear requiring a RG analysis to determine the 
scaling behavior in the neighborhood of the depinning 
transition. This has been done using different approaches 
by Natterman et al |48) and Narayan and Fisher |^9],[50| . 
The order parameter of the depinning transition is the 
average interface velocity. For F < F c we have v = 
while for F > F c the interface moves with the average 
interface velocity v ~ (F — F C )P , with 



/3 =!/(*-£), 



(5) 



where v is the correlation length exponent. 

The upper critical dimension is d c = 4 and the scal- 
ing exponents can be obtained through a functional RG 
analysis. Below the upper critical dimension it results 
that 



c 



A-d 



(6) 



Moreover, near the critical state the static susceptibil- 
ity and crorrelation length diverges according to \ ~ 



2 



(F — F c ) 7 and £ ~ (F — F c ) respectively, where the 
scaling exponents 7 and v satisfy the scaling relation 



Alternative to the continuum equation ([j]) one may 
consider the following discrete model J|^] . In a d + 1 hy- 
percubic discrete lattice each site is labeled by index i 
and height hi. Discrete time steps t — 0, 1, . . . are taken. 
The force acting on each site and the evolution of the 
interface height hi(t) are given by 

F i = X] h i ~ 2dhl + F + Viihi), 

II II 

hi{t + l)-hi(t) = G(Fi), (8) 

where <d(x) is the Heaviside unit step function, nn de- 
notes that the sum runs over the 2d nearest neighbors 
and J2nn hj — 2dhi is a discretized Laplacian. On each 
step t and at every site where Fi > (O(i^) = 1) the 
interface is advanced in parallel by one unit. Thus, only 
at sites where Fi > the interface is active. The average 
interface velocity is thus given by (in non-dimensional 
variables) 

W = Pa, (9) 

where p a is the density of active sites. In the continuous 
representation this expression is equivalent to 

(d t h{x,t)) = ( Pa (x,t)}, (10) 

where p a (x, t) is the coarse-grained density of active sites. 

LID models are Abelian, in the sense that the order 
in which sites advance is not important p2| . If a site i 
is active then it will transfer energy to its nearest neigh- 
bors, which at the same time may become active, and so 
on, an avalanche is generated. It is thus possible that at 
certain time step t there will be more that one active site. 
These sites will be updated in parallel according to the 
evolution rules described above. However, the process 
of toppling can never transform any active site, different 
from itself, in inactive and, therefore, the other active 
sites will remain active. On the other hand, the energy 
transferred to its neighbors is constant, independent of 
the total force at this site. Hence, the order in which 
these sites are updated is not important. 

III. CONSTANT RATE LID MODEL 



linearly with time. Hence, if we are trying to map any 
of these systems into a LID model we can not assume a 
constant force, which will not be in correspondence with 
the picture described above. 

Motivated by this fact we propose a LID model where 
the force increases linearly with time. However, in order 
to reach an stationary state we most include a restor- 
ing force which balance the external field. If the force 
increases linearly with time (F = ct) then after certain 
time it will overcome the critical force F c and the inter- 
face will start moving with a velocity v(F — ct), which 
is time dependent. Hence the system could not reach a 
constant velocity stationary state. This problem can be 
solved adding a restoring force which balance the exter- 
nal driving and leads to an stationary state. This can be 
done in different ways. For instance, one may consider a 
local restoring force linear in h resulting 

F = ct-eh, (11) 



In this case the interface is driven by an external force 
which increases at rate c and a restoring force which pull 
the interface to the substrate. The coefficient e is a mea- 
sure of the strength of the restoring force. Now, suppose 
that at the initial state h(x, 0) = 0. With increasing time 
the external force F(t) — ct will increase. At certain time 
F > F c and the interface will start moving with an av- 
erage velocity v, i.e. (h) — vt. On the other hand, just 
when h is finite the restoring force will start pulling the 
interface to the substrate h — 0. On average the restoring 
force will be Fn(t) = —evt. An stationary state will be 
obtained when the restoring and pinning forces balance 
the external driving force (ct — F c + evt) resulting, in the 
stationary state, the average interface velocity v = c/e. 
We can also have a model with e = but with the inter- 
face pinned at the boundary, i.e. ft, = at the boundary. 
In this case it is expected that the interface develop a 
parabolic profile so that the surface tension balance the 
action of the driving force. 

Alternative, instead of the local restoring force in eq. 
( |TT| ) one may consider a golbal restoring force as follows 

F = ct-eJ -^h(x',t), (12) 



In most systems which are expected to exhibit SOC 
the external field, instead of being constant, increases 
linearly with time. For instance, in sandpiles grains are 
usually added at constant rate so that the total number 
of grains added to the pile up to time t increases linearly 
with t. A more evident example is found in magnetic 
noise measurements in ferromagnetic and superconduc- 
tor materials where the applied magnetic field increases 



The qualitative analysis developed below will be also 
valid for this case. Moreover, we have recently shown 
that both restoring forces lead to the same critical be- 
havior |5l| ]. In the following we will only consider the 
local restoring force. 

Next we proceed to show how different SOC systems 
can be mapped into a constant rate LID model. 
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A. The Barkhausen effect 

An inmediate realization of the equation of motion (|l|) 
with a local force given either by eq. ( p"l| ) or ([l2]) is the 
dynamics o domain walls. It is well known that domain 
walls in ferromagnets move following an irregular motion 
in response to changes in an applied external magnetic 
field, leading to discrete jumps in the magnetization, a 
phenomenon known as the Barkhausen effect. In fact 
Urbach et al H considered a model where the domain 
wall dynamics is described by a LID model, with either 
a local o global restoring force. More recently Cizeau et 
al H have shown that the model by Urbach et al is only 
valid if long-range dipolar interactions are neglected, an 
approximation which is valid in soft magnetic materials 
J7| . Hence, the dynamics of a domain wall in a soft mag- 
netic material can be described by a LID model driven 
at constant rate. 



B. Sandpile models 

In cellular automaton sandpile models a discrete or 
continuous variable Zj, height or energy, is defined in a 
rf-dimensional lattice. The dynamical evolution of Zj is 
defined by two evolution rules: adding and toppling. Dif- 
ferent adding and toppling rules may be defined leading 
to different models E]P]. In particular we consider the 
following rules 

• adding: on each step each site receives a grain from 
the driving field with probability c; 

• toppling: if at certain site Zj > z c = 2d — 1 then 
zi — > Zi — 2d — e + rji and z nn — » z nn + 1 at the 2d 



nearest neighbors; 

e is the average dissipation rate per toppling event and r\i 
is a noise, which may have different origins. If e > the 
toppling rule is non-conservative and the model have bulk 
dissipation. On the contrary, if e — the toppling rule is 
conservative and one has to assume open boundary con- 
ditions to balance the input of grains from the external 
field. For instance, Chessa et al Q considered a non- 
conservative sandpile model where the toppling site loses 
its energy with probability p without transferring it to 
its neighbors. This corresponds to an average dissipation 
rate per toppling event e = 2dp, while r)i will reflect the 
stochastic nature of the local dissipation. Another exam- 
ple is found in random threshold models |52] . In this case 
after each toppling event a new random threshold is as- 
signed (a new critical slope in the model by Christensen 
et al [p2|), which is equivalent to introduce the noise rji. 
On the other hand, at a coarse-grained level the noise r)i 
will appear due to elimination of microscopic degrees of 
freedom. Hence, at this level, it will be present no mat- 
ter if the model is stochastic or deterministic. In other 
words, we expect that deterministic and stochastic mod- 
els belongs to the same universality class. 



In the original BTW model e = and dissipation takes 
place at the boundaries, while grains are added only after 
all sites become stable, i.e. z 7 ; < z n at all sites. However, 
it has been shown jjj],[19 '. 
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that the BTW model and 
the model with bulk dissipation lead to the same critical 
behavior provided e ~ L~ 2 , where L is the lattice size, 
and c — > 0. In the limit c = + we have separation of time 
scales between avalanche duration and energy addition, 
as assumed in the BTW model. Moreover, as discussed 
above, the noise rji will appear at a coarse-grained level. 



Paczuski and Boettcher [ 
critical slope sandpile mode 



24}] noted that a particular 
can be mapped into a LID 



model where hi(t) is the number of toppling events at 
site i up to time t. Their analysis was limited to one 
dimension but can be extended to larger dimensions and 
models with critical height rules [^,g6|. If at t = we 
have Zi = at all sites then, at time step t, Zi and hi are 
related via 



^ hj - 2dhi +ct- ehi + Vii^i), 



(13) 



where J2nn hj gives the number of grains received from 
nearest neighbors, 2dhi the number of grains transferred 
to nearest neighbors, ct the number of grains received 
from the driving field, and —ehi the number of dissipated 
grains. When a site topples hi — > hi + 1 and therefore the 
interface profile hi(t) will always advance in the positive 
direction. Since a site topples only when Zi > z c then 



hi(t+l) - hi(t) = 0(z, - z c ), 



(14) 



Instead of follows the evolution of the variables Zj, 
which are equivalent to the force in LID problems, one 
may follow the evolution of hi, which is equivalent to 
the interface height in LID models. Notice that the evo- 
lution equations ([l3|) and (|l4|) for sandpile models are 
equivalent to the discrete variant of LID models in eq. 
|8|), taking Fi = z» — z c . Hence, after coarse-graining, we 
obtain an equation of motion for h(x,t) like (0) with F 
given by eq. ([ll]). Actually there will be an additional 
term — z c in the right hand side of eq. (0) but it does not 
carry any changes in the critical behavior so that one can 
work without considering this term. 



IV. PHASE DIAGRAM AND SCALING 
EXPONENTS 

We have shown that different SOC models can be 
mapped into a LID model where the force increases at 
constant rate and a restoring force pulls the interface to 
the substrate. In this section we solve this model us- 
ing previous results for the constant force variant. More 
precisely, we investigate the dynamics of an interface de- 
scribed by the equation of motion ([!]) with F given by 
eq. ©, i.e. 



\d t h = TV 2 h + ct-eh + rj(x, h). 



(15) 
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As we discussed above, in this case the interface will never 
be pinned but moves with a finite average velocity v. It 
is thus expected that the dynamics will be similar to 
that observed in constant force models above the critical 
force. Based on this supposition we will perform a suit- 
able change of reference in order to obtain an equation 
similar to that for the fluctuations around the average in 
constant force models, eq. (0). 



A. case e = 

Let us first consider the case of boundary pinning 
which corresponds with e = and h = at the boundary. 
We look for a solution in the form 



h(x,t) = h (x,t) +y(x,t), 



(16) 



so that when we substitute this expression in eq. (|l5| ) 
the constant rate force will be replaced by a constant 
force. This constraint will be satisfied taking ho(x,t) as 
the solution of the problem 



V 2 h Q + ct = F. 



(17) 



with the boundary condition ho = 0. The solution of this 
equation can be easily obtained assuming radial symme- 
try ho = ho (r, t) with < r < R, where R is the system 
radius. It results that 



h (r, t) 



r 2 

l -R> 



ct-F 
2Td 



-Rr 



(18) 



We have also analyzed the case of free boundary condi- 
tions dho(R,t)/dr = 0, resulting that eq. (|l7| ) have no 
solution. In this case the system will not reach an sta- 
tionary state. This result supports our previous analysis 
about the need of a restoring force which balance the 
external driving of the interface. 

Substituting eq. ([16]) in eq. ( |l5[ ) and taking the limits 
r <C R and ct 3> F we obtain eq. (0) with 



Oho 
dt 



cR? 
2Td' 



(19) 



Notice that assuming r <C R and ct ^> F we have ap- 
proximated h by vt + y in the quenched noise rj(x,h). 
However, the information contained in ho is not com- 
pletely lost because we have the constant force F in the 
right hand side of eq. (^J). These approximations will 
be valid for very long times (stationary state) and far 
from the boundary, v is then the average interface ve- 
locity in the stationary state, which depends on system 
size R. When increasing system size the bulk increases 
faster than the surface. Hence, since the interface is only 
fixed at the boundary the concentration of points where 
the interface is fixed will decrease with increasing system 
size resulting on an increase of the interface velocity. 

The constant force F has not been specified. It 
will be obtained self-consistently imposing the constraint 



(y(x,t)) — 0. In constant force models F is a fixed 
parameter while v is obtained self-consistently from the 
equation of motion. However in the present model v is 
given by eq. (^|) while F is the undetermined param- 
eter. From the constant force variant it is known that 
a force F (F > F c ) gives an average interface velocity 
v ~ (F — Fc) 13 . Hence, in the constant rate model, to 
obtain a velocity v we should have the force 



F(v) = F c + const, v 



1//3 



(20) 



In this way, F(v) is obtained self-consistently imposing 
that the interface moves with an average velocity given 
by eq. (|l9|) . In spite of this difference both models have 
the same critical behavior. 

To reach the critical force F c we must fine-tune v to 
zero. For v > there is a characteristic length £ ~ v~ v l^ '. 
Now, the average interface velocity scales as v ~ cR 2 and 
therefore when we take the thermodynamic limit R —> oo 
we must fine-tune c to zero to reach the critical state. If 
£ 3> R then the system size R will be the only character- 
istic length and the system will be in a critical state. This 
condition is satisfied if c = 0(R~ 2+V '@) when R — > oo. 

In the case of boundary pinning the average density of 
active sites, computed from eqs. (|l0|), ( |l6| ) and (|l8|), is 
given by 



(Pa(r,t)) 



r 2 \ c 
#2 ) 2Td 



R z 



(21) 



This expression is identical, except for some constant fac- 
tors, to that obtained in the field theory by Vespignani 
et al It just reflects the balance between the driv- 
ing force and the boundary pinning (the driving field 
and boundary dissipation in sandpile models). As one 
can see from eq. (|2l|) the average density develops a 
parabolic profile, which is a consequence of the pinning 
at the boundary. This parabolic profile has been observed 
in recent numerical simulations by Barrat et al |53| . 



B. Case e > 

Now we consider the case of a local restoring force. We 
again look for a solution in the form of eq. ( |l6| ) which 
the same constraints, i.e. it replaces the constant rate 
force by a constant force F. This constraint is satisfied 
if ho(x,t) is the solution of the problem 



V 2 h - eho + ct = F. 



(22) 



with the corresponding boundary conditions. Indepen- 
dently of the boundary conditions assumed, the solution 
of this equation exists. If one assumes the fixed boundary 
condition then one will obtain a parabolic profile as in the 
e = case. On the contrary, we will consider periodic or 
free boundary conditions. In this case we obtain 
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h (x,t) = ——. (23) 

Substituting eq. (^3|) in eq. (|l^) and taking the limit 
vt 3> F/e we obtain 

Xd t y = TV 2 y - ey + F - \v + r)(x, vt + y). (24) 

with 

v = -^. (25) 

Again assuming vt ho(x) we have approximated h by 
vt + y in the the quenched noise rj(x,h). Eq. (|24|) is quite 
similar to that obtained for the fluctuations around the 
flat co-moving interface in constant force LID models, 
eq. ([!]). The only difference is found in the term —eh, 
which accounts for the local restoring force with strength 
e. However, in the limit e — * the system will show the 
same critical behavior as in constant force models. The 
term — ey only affects the linear part of the bare prop- 
agator while non-linear terms, which are responsible for 
loop-corrections, remain identical. The critical exponents 
z and £ obtained from the RG analysis will thus be the 
same as those obtained for the constant force case. How- 
ever the phase diagram will show a complex structure. 

For v > and e > the system is driving out the crit- 
ical state. Now in addition to the characteristic length 
£ ~ v^lV we have another characteristic length associ- 
ated with the restoring force £ e ~ e - "", where v e is cal- 
culated below. The phase diagram (v,e) will have dif- 
ferent regions depending on the ratio between these two 
characteristic lengths. We then define the characteristic 
velocity v e as the velocity where these two characteristic 
lengths become identical. Thus taking £ ~ £ c we obtain 
that v e ~ e^ 8 with 

& = ^(*-0- (26) 

In the region v > v e (£ > £ c ) £ is the only character- 
istic length, as in the case e = 0. However, we cannot 
say that the effect of the restoring force disappears be- 
cause the average interface velocity depends on e (see eq. 
(p5|)). There will be a complete equivalence if we take 
e ~ R~ 2 . In other words, a model with bulk dissipation 
with v > v c and e ~ R~ 2 is equivalent to a model without 
bulk dissipation and open boundaries. This equivalence 
was already pointed out by Vespignani et al using MF 
jjJI and field ]l9) theory approaches. 

On the contrary, in the region v < v e (£ < £ e ) £ e is the 
relevant characteristic length making the difference with 
the e = case. In this region the static susceptibility, 
which characterizes the linear response of the system, is 
given by 

X (fc,e) = l/(r off fc 2 + e). (27) 

where T s g includes loop corrections to T. At the critical 
state e — > we have 0) ~ fc~ 2 as in the constant force 



case. On the other hand, for k — > we obtain \ — V e so 
that x ~ e_7e with 

% = I- (28) 

Eq. (|27j) can be thus written as 

X (k,e)~k- 2 f(k&, (29) 

where £ c ~ e~"' . For k^ e large we should obtain x(fc, e) ~ 
fc~ 2 so that /(a;) ~ 1 for large cc. On the other hand, for 
kt; small we should now obtain x(fc, e) ~ e _7e and there- 
fore f(x) ~ x 7 '/"' for x small with 

^ = 2. (30) 

From eqs. ( P8|) and (|3^) we thus obtain 

= 1/2. (31) 

The exponents j e and v t results different to the exponents 
7 and v but their ratio is the same, as one can see from 
eqs. (^) and ([30]). This will carry as a consequence that 
some exponents, the avalanche exponents for instance, 
measured in the region v < v c can be extrapolated to the 
region v > v c . 

Finally, in the case of a local restoring force the average 
density of active sites is given by 

(p (r,i)> = ~. (32) 

Again this result is identical to the one obtained by 
Vespignani et al [|l9| using a field theory for sandpile mod- 
els. It reflects the balance between the driving and the 
local restoring forces (between driving and local dissipa- 
tion in sandpile models). Moreover this flat profile has 
been observed by Barrat et al |53[ in numerical simula- 
tions of a sandpile model with local dissipation. 

The constant rate LID model thus describe a wide vari- 
ety of SOC systems, like the Barkhausen effect and sand- 
pile models. It reproduces previous field theory predic- 
tions for sandpile models, which we now know are also 
valid for other models. Hence, the constant rate LID 
model provides a unifying point of view of SOC phenom- 
ena. 



C. Avalanche exponents 

In the preceding section we have shown that LID mod- 
els driven either at constant force or rate belong to the 
same universality class. In particular they share the same 
roughness and dynamic scaling exponents, which can be 
estimated using previous RG calculations for the con- 
stant force LID model. Now we proceed to show how 
other scaling exponents can be obtained using £ and z. 
For instance, we are going to compute the avalanche ex- 
ponents, which are often measured in numerical simula- 
tions of sandpile models. 



G 



Let s be the avalanche size and T its duration, which 
are distributed according to P(s) and P(T), respectively. 
Just at the critical state one expect that these distribu- 
tions satisfy the power law behavior P(s) ~ s~ Ts and 
P(T) ~ T~ Tt , where r s and r t are the avalanche distri- 
bution exponents, reflecting the unexistence of character- 
istic values for size and duration of the avalanches. How- 
ever, for finite e and v characteristic cutoffs of avalanche 
size s c and duration T c will appear. These cutoffs will 
scale with the correlation length (either £ or £ e depend- 
ing on the model and on the region of the phase diagram) 
as s c ~ £ D and T c ~ £ z , where D is the avalanche frac- 
tal dimension. Near the critical state the distributions of 
avalanche size and duration will thus satisfy the scaling 
laws 

P(s) ~ «- T '/i(*/*c). ~ T-^ 9l (T/T c ), (33) 

where /(x) and g{x) are some cutoff functions with the 
asymptotic behaviors /i(x), g\{x) ~ 1 for i C 1 and 
/ x (x),t/i(x) < 1 for x > 1. 

The exponents r s , rt, D and z are not all independent. 
Since s ~ T Z ' D then the condition / dsP(s) = J dTP{T) 
implies 



(t s - \)D = { n - i)z. 



(34) 



Another scaling relation can be obtained taking into 
account that the short-wavelength static susceptibility 
scales as the mean avalanche size Jl2| . From eq. ([33]) we 
obtain (s) ~ £_( 2 -^) D while X ~ £ 2 . Thus, setting 

X ~ (s) we obtain 



(2 - t,)D = 2. 
Then from eqs. (133) and (55|) it results that 



r, = 2 , 

D 



T t = 1 



D 



(35) 



(36) 



Finally there is a scaling relation which relates the 
avalanche dimension exponent D with the roughness 
exponent £. Below the upper critical dimension the 
avalanches are compact objects ||^] and therefore s ~ 
Ahr d , where Ah is the characteristic fluctuation of the 
interface during the avalanche and r its characteristic 
linear extent in the d-dimensional substrate. Then, since 



Ah ~ and s 



one obtains fB2| 



D = d + C 



(37) 



Above the upper critical dimension the avalanches are 
no more compact and D = d c = 4 |32|, resulting the 
mean-field exponents r s = 3/2 and r t = 2. 

These scaling relations are independent on the way the 
interface is driven. They have been obtained assuming 
that x ^ (s), s ~ T D ' Z , there is a characteristic length 
£ and avalanches are compact objects (in the interface 
dynamics representation) . 



V. DISCUSSION 

In this section we proceed to compare our results with 
experiments, numerical simulations and previous theo- 
retical works in the literature. 

The average density of active sites in the stationary 
state was found to be identical to the one obtained in the 
field theory by Vespignani et al |is[| , which in principle 
was developed for sandpile models. Moreover, they pre- 
dicted an upper critical dimension d c = 4 in agreement 
with our result. Our analysis thus reveals that their ap- 
proach is also valid for other SOC models, for instance 
for LID models driven at constant rate. The advantage 
of our approach is that it makes easier the use of momen- 
tum space RG analysis to obtain the scaling behavior in 
the neighborhood of the critical state. Using previous 
RG calculations for constant force LID models we com- 
pute the complete set of scaling exponents, including the 
avalanche exponents often measured in numerical simu- 
lations. 

The upper critical dimension d c — 4 is also consistent 
with previous RG analysis by Dfaz-Guilera |T(| of the 
corresponding Langevin equations for the height of sand 
columns Zi (instead of hi) in the BTW model. Their re- 
sults are also consistent with a dynamic scaling exponent 



d + 2 



(38) 



an expression which was previously suggested by Zhang 
p5[ . This dynamic scaling is however in contradiction 
with the one obtained here (see eq. @. In his for- 
mulation Dfaz-Guilera considered a coarse graining-noise 
r/(x, t). In the fast time scale (that of the evolution of the 
avalanches) they assumed a columnar noise uncorrelated 
in space, i.e. 



(rj' (x, t)rj'(x' , t') ~ 5 d (x — x'). 



(39) 



On the contrary, in our approach rj{x, h) is a quenched 
noise, a new random variable is selected only when the 
interface at site x advances, i.e. when the site x is active. 
Our choice is more appropriate because it reflects the fact 
that the noise field is frozen in regions where there is no 
dynamical evolution, similar to the multiplicative noise 
assumed in absorbing-state phase transitions Jlgf ]. 

The estimate of the avalanche scaling exponents of 
sandpile models has been a very difficult task. In two 
dimensions we count with a real space RG approach 
by Pietronero et al flq| , which predicts that undirected 
Abelian sandpile models, either deterministic or stochas- 
tic, belong to the same universality class. They also 
provided an estimate of the avalanche size exponent 
T,,. = 1.253. On the other hand, Priezzhev et al |14| have 
obtained t s = 6/5 = 1.2 in two dimensions, analyzing the 
structure of avalanches in the undirected Abelian sand- 
pile model. Since the BTW and Manna model are both 
Abelian |13| this result will be valid for both models, 



7 



which are thus expected to belong to the same universal- 
ity class. Moreover, this exponent is identical to the one 
obtained by De Menech et al [|o| for the BTW model, 
using multifractal rather than finite size scaling. These 
estimates can be compared with the one obtained within 
our approach. In two dimensions we obtain r s = | = 1.25 
(see table [j]) . This value is in very good agreement with 
the theoretical estimates reported above. The advantage 
of our approach is that we go beyond two dimensions 
and compute the scaling exponents in three dimensions. 
Moreover we have determined the upper critical dimen- 
sion, d c — 4 above which MF exponents are correct. 

Numerical simulations of undirected sandpile models 
are generally performed with dissipation only at the 
boundary and assuming time scale separation. In this 
case the system self-organizes into a critical state where 
the system size is the only characteristic length, i.e. 
£ ~ L. Moreover, it is observed that the susceptibility 
scales as \ ~ L 2 . This scale dependency has been demon- 
strated for Abelian models in one dimension |13| using 
the Abelian symmetry of the operators algebra. It is ex- 
pected to be valid in any dimension as a consequence of 
conservation |l^,|l9|, understanding conservation as the 
balance between the input and output of grains from the 
system. This scaling dependence has been obtained here 
using the mapping of sandpile models to the constant 
rate LID model, and it is a consequence of the balance 
between the driving field and the restoring force. 

The scaling exponents obtained from numerical sim- 
ulations of constant force LID models, the Manna 
model (the prototype of stochastic sandpile model) and 
Barkhausen effect in soft-magnetic materials are shown in 
table [I]. Some experimental estimates for the Barkhausen 
effect measure in soft magnetic materials, where dipolar 
interactions can be neglected |J7[. are also displayed. Our 
estimates using eqs. ( J36| ) and (|37]) are also shown for com- 
parison. The agreement is very good, specially in three 
dimensions where the RG calculations are expected to be 
more precise. Moreover, we can also observe that expo- 
nents obtained for the three different SOC models in two 
dimensions are very close, suggesting that they all belong 
to the same universality class. In this analysis we have 
not included numerical estimates for the avalanche expo- 
nents of the BTW model. Numerical simulations of this 
model ^7|^8| lead to contradictory results. This discrep- 
ancy may be associated with the existence of multifractal 
rather than finite size scaling p0| . 

This work does not cover the whole classes of SOC 
models. We have excluded, for instance, directed Abelian 
sandpile models, critical slope models |57| an d non-linear 
interface depinning models ]42|] . In the case of directed 
Abelian sandpile models the critical exponents have been 
obtained using the Abelian symmetry of the algebra of 
operators jL3|. In the case of critical slope models, in- 
spired on the work by Paczuski and Boettcher p4j , one 
can perform a similar mapping into a LID model fl57| . 
Finally non-linear interface depinning models belongs to 
a different universality class described by the KPZ equa- 



tion 



VI. SUMMARY AND CONCLUSIONS 

We have shown that undirected Abelian sandpile mod- 
els and the domain wall motion in magnetic materi- 
als where dipolar interactions are negligible can all be 
mapped into the constant rate LID model and therefore 
display the same universal behavior. The constant rate 
LID model was found to exhibit the same critical behav- 
ior as constant force LID models above the depinning 
transition. The scaling exponents were thus obtained us- 
ing previous RG analysis for constant force LID models. 

In this way we have shown the existence of universal 
behavior in a vast class of SOC models. We have pro- 
posed the LID driven at constant rate as the prototype of 
this universality class, because it allows the use of contin- 
uous approaches which are more treatable by analytical 
tools than the discrete analysis required in most cellular 
automaton SOC models. 
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d Model 



Tf 



D 



Ref. 



1 LID 
RG 

2 LID 
Manna 
Manna 
BHN 
BHE 
RG 

3 LID 

Manna 
RG 



1.29(2) 

1.273 

1.27(1) 

1.26(4) 

1.28(2) 

f = 1.25 



1.40 



1.50(1) 
1.40(5) 
1.5(1) 
±2 « 1.43 



1.42(3) 
| « 1.33 

1.58(4) 

1.500 

1.50(2) 



£ « 1.56 



1.75 



1.75 

16 ~ 
9 ~ 



1.78 



2.25(1) 
2 

2.75(2) 

2.750 

2.73(2) 



| w 2.67 

3.34(1) 
3.33 

f » 3.33 



35 



35 
27 
54 

r 



a 



27 



TABLE I. Scaling exponents for constant force LID models 
(LID), the Manna d-state model (Manna), and those obtained 
from experiments (BHE) and numerical simulations (BHN) of 
the Barkhausen effect. Results obtained here using RG esti- 
mates are shown for comparison. 
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